arXiv:1509.02790vl [math.NA] 9 Sep 2015 


1 


On the Hierarchical Preconditioning of the Combined 

Field Integral Equation 

Simon B. Adrian^Francesco R Andriulli^, and Thomas F. Eibert^ 

^Lehrstuhl fiir Hochfrequenztechnik, Technische Universitat Miinchen, Munich, Germany 
^Institut Mines-Telecom / Telecom Bretagne, Technopole Brest-Iroise, Brest, France 


Abstract —This paper analyzes how hierarchical bases precondi¬ 
tioners constructed for the Electric Field Integral Equation (EFIE) 
can be effectively applied to the Combined Field Integral Equation 
(CFIE). For the case where no hierarchical solenoidal basis is 
available (e.g., on unstructured meshes), a new scheme is proposed: 
the CFIE is implicitly preconditioned on the solenoidal Helmholtz 
subspace by using a Helmholtz projector, while a hierarchical non- 
solenoidal basis is used for the non-solenoidal Helmholtz subspace. 
This results in a well-conditioned system. Numerical results corrob¬ 
orate the presented theory. 

Index Terms —integral equations, preconditioning 

1. Introduction 

Scattering and radiation problems in electromagnetics are com¬ 
monly solved by leveraging the Electric Field Integral Equation 
(EFIE) in the framework of boundary element methods. The 
EFIE, however, is not well-conditioned: whenever the frequency 
decreases or the discretization density increases, the condition 
number grows [1], [2]; these effects are often called low- 
frequency and dense-discretization breakdown. In addition, the 
EFIE is not uniquely solvable at resonance frequencies [3]. 

The low-frequency breakdown has been solved in the past by 
using a quasi-Helmholtz decomposition of the current, such as the 
loop-star/loop-tree decomposition, allowing for a separate rescal¬ 
ing of the vector and the scalar potential in frequency [2]. The 
dense-discretization breakdown was remedied by the introduction 
of hierarchical basis preconditioners [4]-[7]. Since hierarchical 
basis preconditioners are usually based on a quasi-Helmholtz 
decomposition, they cure the low-frequency breakdown as well. 
The interior resonances have been removed by solving a linear 
combination of the EFIE and of the Magnetic Field Integral 
Equation (MFIE): the Combined Field Integral Equation (CFIE) 
[3]. Since, however, the CFIE is a linear combination of the EFIE 
and of the MFIE, it is not free from the low-frequency and the 
dense-discretization breakdown. For this reason, efforts have been 
made to extend the hierarchical basis preconditioners to the CFIE 
[8]. 

The conference contribution [8] reported on the fact that 
the application of a hierarchical loop (h-loop)/hierarchical non- 
solenoidal (h-n-sol) basis preconditioner to the CFIE resulted in 
a well-conditioned equation. Given that the h-loop preconditioner 
is not capable of preconditioning the EFIE [7], [9], we found the 
result surprising and decided to perform a theoretical investigation 
of the problem. This lead to the current letter whose novelty 
content is twofold: (i) we show on a theoretically sound basis 
that a direct application of the hierarchical basis preconditioner 


to the CFIE is possible when the h-loop [4] or the 3-point 
hierarchical loop (3ph-loop) functions [7] are used as solenoidal 
basis, (ii) We propose a new hierarchical left preconditioning 
scheme that differently from any other hierarchical approach to 
the preconditioning of CFIE (including those communicated in 
[8] and detailed in [10]) can provide optimal preconditioning 
both on structured and unstructured meshes. Numerical results 
will show the effectiveness of the new scheme for canonical 
and realistic examples. Very preliminary results have appeared 
as conference contribution [11]. 

II. Notation and Background 

We consider an electromagnetic wave W) impinging on a 
Perfectly Electrically Conducting (PEC) scatterer Q with bound¬ 
ary F = dQ, which is embedded in a homogeneous medium with 
permittivity s and permeability jd. The impinging wave excites an 
electric surface current density j generating the scattered wave 
(E\W). Thus the total field is (E,H) = (E^aE\WaW) 
and on F the tangential component of the electric field satisfies 
the boundary condition h X E = 0 . 

The surface current density j is recovered from the tangential 
component of the incident electric field by solving the EFIE 

j = -hxE^ and from the tangential component of the magnetic 
field by solving the MFIE A\j = {112 A *K) j = rfh X W 
on the surface F with wave impedance r] = where we 

assumed that the current is normalized such that j = rjh X H; 
the operator T' = 7a + To is the EFIE operator with the 
vector potential TaJ = hxik scalar potential 

nj = -hx^Vr ^^yr-j(r')dr', and M = IH+K is 
the MFIE operator with ’Kj = —hx{Vr X 
where k is the wavenumber and h is the exterior surface normal 
vector. The CFIE is given by Cj = —ahxhxE^+ri{\—a)hxW, 
where C = ah xT A {I - a)A\ is the CFIE operator with 
a €]0,1[. 

Following a Galerkin approach, the CFIE is commonly dis¬ 
cretized using Rao-Wilton-Glisson (RWG) functions f i € Xrwg 
as source and testing functions. Hence, the CFIE system reads 

Zci = vc (1) 

where = ceZ^ + (1 - a)T]ZM, j ~ [']jfj’ 

Vq = av£ -l- (I — a^Tj \/]y[ with = Za + ^o? 

[ve]. = - ifi, nxhx E^)j ^2 and [Za] = h x 

[^Aij = [Zm].. = 

[vy,]. = {fi,hxW)^2. 
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It has been pointed out that testing the MFIE operator with 
standard RWG functions is not consistent with the correct 
discretization framework of the operator and that one way of 
obtaining a correct and conforming dis^etization would be to 
use Buffa-Christiansen (BC) functions f i as testing functions 
[12]. In this way,^he conformingly discretized MFIE system reads 

[^Aii = ^/y)i :,2 and [i/m], = (n Xn X 

Accordingly, the modified CFIE matrix reads 


\Zc = aZ-E + (1 


■ jZm, 


( 2 ) 


and Vq = av^ + {I - cr)7/GffG^^~ ^\/m with the Gram matrices 

[Gff].. = and = [n xf f j) 

It is well known that the EEIE is ill-conditioned with 
cond(ZE) = O [ij where h is the average edge length 

of the mesh. Hierarchical basis preconditioners that cure the 
EFIE ill-conditioning consist of two sets of basis functions, 
hierarchical solenoidal (h-sol) and h-n-sol functions forming a 
quasi-Helmholtz decomposition of ^Crwg- 

Before we study the hierarchical preconditioners, we consider 
the loop-star quasi-Helmholtz decomposition of the EFIE. Let 
Ai be the loop functions and 27/ be the star functions. Then we 
have j ^ ['a] .^7 + SS [k]j Zj. Since 

the loop and star functions are denned as linear combinations 
of RWG functions, we can define transformation matrices A and 
X such that / = Ai\ X/^. We define T = |^/l/V^ XV^j 
and the EFIE system matrix in the loop-star basis is given by 

Zeqh •= . In the low-frequency limit, 

^E,qH becomes block diagonal, that is, Z\y and Z^a vanish; 
thus, Zaa and Z^e can be identified as the sources of the 
ill-conditioning. In fact, it has been shown that the matrices 
Z'aa,o = limit ^0 Zaa and Zee,o = hmyt^o Z^y can be expressed 
as [Zaa,o].. = and [Zss,o].. = (cr,-, 

where T/ are the nodal functions (i.e., Ai(r) = 1 for r e Ui, 
Ai(r) = 0 for r e nj ^ Ui, where nt c L is the iih node of the 
mesh), (Ti = divX/, and Ai = V X hAi. Moreover the operators 
are defined as AVf = hr ' V X ^ ^r'f{r')dS{r') 

and = fr W=^/('-')d5(r') [7], 


III. Spectral Analysis and Preconditioning 

The operators AV and A^ are known to suffer from the dense- 
discretization breakdown [7]. If we precondition these operators 
successfully, we can precondition the EFIE. Different hierarchical 
bases have been introduced in the past. For preconditioning 
To, the bases developed are all equivalent, and we choose 
as a representative the basis developed in [5]. For 7 a, two 
hierarchical bases have been presented: the h-loops Aui [4] 
and the 3ph-loops At/ [V]. Only the latter can cure, however, 
the dense-discretization breakdown of the EFIE [7]. Both Ah^ 
and At/ require a structured mesh: a mesh that is obtained 
by iteratively refining, for example dyadically, J times a coarse 
mesh. Accordingly, Ah^ and At/ are defined on J hierarchical 
levels (for details, see [7]). These functions are related to scalar 
functions by Ani = V X hAui and At^ = V X hAji, where 
Th/ are the hierarchical nodal (h-nodal) functions presented by 
Yserentant [13] and Tt/ are the 3-point hierarchical nodal (3ph- 
nodal) functions by Stevenson [14]. As we did for the loop basis. 


we can define transformation matrices Ah and At for the h-loop 
and 3ph-loop functions, respectively. 

The operator AV maps from the Sobolev space 
to and induces an inner product on (for 

a definition of the Sobolev spaces, see [15]). The 3ph-nodal 
functions Aji are /^-stable when rescaled appropriately, and 
therefore, the 3ph-loop functions can precondition the vector 
potential [7] resulting in cond (^D^^AjZaAjD^j^^ < 1, where 
is a diagonal matrix with entries [^a] - “ where /(/), 

i e {0 ...,/}, is the number of the level of the iih 3ph-loop 
function. We note that makes the 3ph-loop functions 
stable, while results in L^-stable and in //^-stable 3ph- 
loop functions. The h-nodal functions, on the other hand, are not 
//^/^-stable [9], and thus, the h-loop functions fail to precondition 
Ta- Likewise, the nodal functions are not //^^^-stable resulting 
in cond (A^ZaA^ < Ijh [7]. We note that an application of a 
Jacobi preconditioner can improve the behavior quantitatively, but 
not qualitatively (see Fig. la). 

Let Xh/ be the h-n-sol functions for which we define the 
transformation matrix Zh. Since A^ is inducing an // inner 
product, we need a rescaling defined by [^e]// “ 

to make divXn E/^“^^^-stable. Then it can be proved that 
cond (DeXj'JZoXhDe) < log^ (see [7] and references 

therein). 

We define the overall transformation matrix 
:= [XD^/V^ XhDeV^], where X can 

be A, Ah, or At. Summarizing, we find for the 
preconditioned EFIE that the loop/h-n-sol basis yields 

^ 1//1 and that the 3ph- 

loops/h-n-sol basis yields cond < 

log^ (l//t^). 

Following these considerations, we study the CFIE. For the 
sake of brevity, the analysis is carried out only for the conform¬ 
ingly discretized CFIE. The findings are, however, the same as for 
the standard CFIE. Because of the identity operator of the MFIE, 
the conditioning of the CFIE is better than that of the EFIE with 
cond (Zc) < Ijh. While the largest singular value of To still 
grows to infinity, the singular values of Ta are shifted from 0 
to 1/2 by X/2. Next, we have to study a Helmholtz-decomposed 
CFIE. Following the argumentation in [16], the compact operator 
TC can be neglected and it suffices to analyze how G/f changes 
the EFIE behavior. 

We notice that when I is discretized with loop, h-loop 
and 3ph-loop functions, the resulting matrix is equivalent to 
the discretization of the Laplace operator in its weak for¬ 
mulation with nodal, h-nodal and 3ph-nodal functions, since 
(Vxn/,-,Vxn /,)^2 = (v/aV/,)^^ [16], The Laplace op- 
erator induces an inner product on 1^, and hence, we need 
-stable basis functions. 

It is well-known that the nodal functions are not -stable; the 
condition number grows with 0(l/h^), and given that the Laplace 
operator Ar is a pseudo-differential operator of order -1-2 and AV 
is a pseudo-differential operator of order -l-l,J;he total order of 
AV -\-Ar is +2, and hence, we have cond (^A^ZqA^ < Ijh^ [16]. 
In other words, loop functions applied to the CFIE result in a 
conditioning worse than when they are applied to the EFIE. The 
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h-nodal functions are not if ^-stable. Yet, for they result in a 
condition number that grows with log^ [13], and thus, we 

have cond ^ llh^. The 3ph-nodal functions 

are if^-stable [14], and thus, we have cond (^D^AjZqAj x 
1. Notice that the matrix D^AjZeAjD^ is ill-conditioned; 
since we use if ^-stable functions, the singular values are now 
accumulating at zero. This does not destroy the well-conditioning, 
as DlA^G«A-rDl is spectrally equivalent to / and thus the 
spectrum is bounded from below. 

As the basis spanned by div Zm has derivative strength when 
rescaled with (since the basis is if“^/^-stable), the basis 
spanned by Zm has integrative strength. Given that I has no 
derivative terms, the matrix D^X^GffXD^ has singular values 
clustering around zero, and since D^X’^ZeX’D^ is spectrally 
equivalent to X, we can conclude that cond ZDs) < 

log^ (l//!^). 

Summarizing, we have for the loop/h-n-sol basis precondi¬ 
tioner cond ^(/f^ ~ h-loop/h- 

n-sol basis preconditioner cond ~ 

log^ (l//z^^, and for the 3ph-loop/h-n-sol basis preconditioner 

^ log2(l//!2). Equivalently 

said, the combination loop/h-n-sol functions does not precondi¬ 
tion the CFIE, while both combinations h-loop/h-n-sol and 3ph- 
loop/h-n-sol basis are a valid preconditioner for the CFIE. 

In practical scenarios, the mesh is typically unstructured and 
thus h-loop and 3ph-loop functions are not available; yet, from 
the presented theory it is clear that the use of loop functions 
is not effective and a different strategy is necessary. Define the 
transformation matrix A := A • Using these orthogo- 

nalized loop functions, we obtain as hierarchically preconditioned 
system 



While this approach is theoretically sound, it lacks practical 
applicability due to the presence of square root matrices. We 
modify the approach by considering a left instead of the split 
preconditioner in (3). Thereby, we obtain 

It has been shown that P\ := A is a projector to the 

solenoidal Helmholtz subspace [17]. Differently from ^ , 

the inverse of A^A can be rapidly computed iteratively: the matrix 
A^A is a graph Laplacian, and black-box preconditioners like 
algebraic multigrid methods can be used to precondition it [16]. 
The left preconditioner is motivated by the fact that when we 
have matrices A, B G and A is symmetric and x^x < 

x^B^ABx < x^x forall x g R'^, we have cond (88^4^ < 1. 

In our case 4 = Zc is not symmetric due to the perturbation by 
the non-static kernel and the operator 7C. Since the perturbation 
is only compact, we don’t expect a significant impact on the 
behavior of the Krylov subspace methods. 


Usually, best results are obtained when the (hierarchical) func¬ 
tions are rescaled by leveraging on a Jacobi preconditioner, that 

is, we use ,. For a far comparison of 

the different bases discussed, solenoidal and non-solenoidal alike, 
each basis is rescaled by using such a Jacobi preconditioner. 

When the conforming CFIE is used, the Gram matrices 
prohibit to efficiently obtain it. For the preconditioner pre¬ 
sented in this work, this problem can be avoided by using 

= '4 ^h^eZh], ., that is, the same procedure as 
for the EFIE can be used. This choice cures the low-frequency 
breakdown (i.e., the condition is independent of the frequency); 
however, a further aligning of the singular values branches 
associated with the solenoidal and the non-solenoidal Helmholtz 
subspace improves the condition number. Tojthis end, by selecting 
a = IIPaZcII^ and p = USDoDoS^Zcjl^, we can define 
M \= P\/a -\- SD^D^S^/P resulting in the system 

MZc/m = Mvq . (5) 

The norms can be estimated rapidly by using a power iteration 
method. 

IV. Numerical Results 

In the following, we compare our new formulation (5) with 
the unpreconditioned CFIE and when standard hierarchical basis 
preconditioners are applied. The standard hierarchical precondi¬ 
tioners we use always consist of the same h-n-sol basis comple¬ 
mented by loop, h-loop, or 3ph-loop functions, respectively. We 
do this for the conformingly discretized CFIE operator Zc, and 
the standard CFIE Zc (denoted as S-CFIE). 

First, we analyzed the dense-discretization stability by using a 
cube with side length 1 m. We used a plane wave excitation and 
fixed the frequency at / = 1 MHz with / = kc/(2jt), where c 
is the speed of light. We varied the average edge length h from 
1.13 m to 0.07 m, and the number of iteration is displayed as a 
function of the spectral index l//i in Fig. la. The results confirm 
the presented theory: in particular, we can conclude from the 
figure that loop functions should not be used with the CFIE, and 
that h-loops can safely be used with the CFIE but not with the 
EFIE. The new formulation (5) performs well both when applied 
to a conformingly discretized CFIE as well as a standard CFIE. 

Next, we verified the frequency stability. Figure lb displays the 
condition number as a function of the frequency /. We see that 
all preconditioned formulations are free from the low-frequency 
breakdown. Figure Ic shows that all preconditioned CFIEs are 
resonance-free. We also see that in terms of the condition number, 
the unpreconditioned CFIE works better than the preconditioned 
counterparts. This is not unexpected, since hierarchical basis 
preconditioners need to be adapted to high-frequency problems. 

To apply our new method to a more realistic structure, we 
employed the model of a space shuttle shown in Fig. 2. The 
electric size of the shuttle is 1/2/1, where A is the wavelength. 
Since the mesh is unstructured, we do not have a h-sol basis, 
and thus, we can complement the h-n-sol basis only with loop 
functions. Table I summarizes our results, where the solver 
tolerance was 1 X 10“^. 
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(a) Dense-discretization stability. 
Fig. 1. Cube: spectral analysis. 


(b) Low-frequency stability. 



(c) Resonance-freeness. 



Fig. 2. Shuttle 


Formulation 

Iterations 

CFIE 

400 

S-CFIE 

417 

CFIE: loops 

253 

S-CFIE: loops 

198 

This work: CFIE 

33 

This work: S-CFIE 

37 


TABLE I 

Shuttle: the number oe iterations eor the dieeerent eormulations with 

SOLVER TOLERANCE 1 X 10“^. 


V. Conclusion 

First, we can conclude that the h-loop functions, which fail 
to precondition the vector potential 7 a of the EFIE, can be 
successfully applied to the CFIE (though the results obtained 
by using 3ph-loop functions are usually better, and they work for 
the EFIE as well). The best results were obtained with our new 
formulation, which differently from all other hierarchical precon¬ 
ditioners, works on both structured and unstructured meshes. 
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